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I. INTRODUCTION 

Pair transport through interacting one-dimensional 
systems sandwiched between two superconductors has 
been the focus of much attention recently. These sys- 
tems not only exhibit interesting physical phenomena, 
such as Andreev reflection and pair transport, but also 
may be incorporated into novel nanoelectronic devices. 
In particular, there have been various proposals for the 
creation and transport of entangled pairs for quantum 
communication. Therefore, it is useful to have well- 
controlled numerical tools for analyzing phenomena as- 
sociated with pair transport through a superconducting- 
Hubbard-superconducting (SHS) system. 

In this paper we present the results of numerical stud- 
ies of a one-dimensional Hubbard chain sandwiched be- 
tween two superconducting contacts. We examine the ef- 
fects of the boundary contacts on the injection and trans- 
port of pairs through the system. We then introduce the 
idea of an extended contact between the superconductors 
and the intervening Hubbard chain which provides for 
improved pair transmission into and out of the Hubbard 
chain. Using this we present a new numerical method for 
determining the effective superfluid weight D of a Hub- 
bard chain. 



II. THE CONTACT MODEL 

There have been various analytic studies of a onc- 
dimensional Luttinger liquid sandwiched between two 
superconductors^*^. Here, we make use of a comprehen- 
sive analysis recently reported by Affleck et. a! 4 . These 



authors integrated out the electron degrees of freedom 
of the superconducting leads, replacing them with effec- 
tive boundary conditions for the Luttinger liquid. In this 
framework, the effective Hamiltonian for an SHS system 
can be written as 

H = H + H 1 (1) 
where Hq corresponds to the Hubbard chain 

L L L 

Hp = -*y] (c\ a c i+ i a + h.c^j + U ^n^mi - ji rij 
i=i i=l i=l 

(2) 

and Hi incorporates the effects of the two superconduct- 
ing leads 



Hi = A L (e^c^c^ + /i.e.) + A R (e l * R 
+ Vi(mi +nn) + V L (n L -\ +n L i). 
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Here, creates an electron of spin up on the £ site. 
The hopping parameter of the Hubbard chain is t, U is 
the onsite interaction energy, and /x is the usual chemical 
potential. As discussed in reference^, the effect of the two 
superconducting leads can be parametrized in terms of 
contact pairing strengths &il,r) and their phases 4>(l,r), 
along with end point scattering potentials V% and Vl ■ 

In the following we will be interested in the symmetric 
case in which Ax, = A# = A and V\ = Vl = V. The first 
term in Hi injects or removes pairs with different phases 
on both ends. In addition, there are effective boundary 
scattering potentials Vnn which arise and play an im- 
portant role in achieving optimal pair transmission across 
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the ends of the Hubbard chain. Integrating out the su- 
perconducting electron degrees-of-freedom can be seen 
as a natural thing to do when the Fermi level lies well 
below the superconducting gaps in the bulk of the super- 
conductors, since the pair fields in the superconductors 
have well-defined average values and negligible fluctua- 
tions. In the Hubbard system, the value of the pair field 
has to be replaced by the fluctuating pair operator. A 
similar approach was used by Kozub^ to study Josephson 
transport through a Hubbard impurity center. 

In their paper, Affleck, el. a! 4 calculate the Josephson 
current and the Andreev reflection probability. For the 
non-interacting half-filled tight-binding chain, they find 
that the maximum transmission probability is 1 (perfect 
Andreev reflection) and it occurs when A = t and V = 0. 
In this case, the Josephson current versus the phase dif- 
ference 4> — 4>r—4>l between the ends has Ishii's sawtooth 
form£. For smaller values of V, the sawtooth is smoothed 
out and starts resembling the Josephson sine shape cor- 
responding to a small Andreev reflection probability (See 
also, Refii). Away from half- filling, Affleck et. al found 
that in order to achieve perfect Andreev reflection, both 
the contact pairing strength A and the boundary scatter- 
ing potential V needed to be tuned to particular values. 
For the non-interacting case, these values are: 

In order to treat the interacting case, these authors em- 
ployed bosonization and renormalization group methods. 
For negative values of U, they showed that the contact 
Hamiltonian renormalizes to the perfect Andreev reflec- 
tion fixed point. Thus, even when the parameters of the 
contact were not fine-tuned for perfect Andreev reflec- 
tion, one recovers the sawtooth form for the Josephson 
current versus the phase difference as the length L of the 
Luttinger liquid increases. However, for positive U, they 
found that the contact Hamiltonian flows away from the 
Andreev fixed point. In this case, as L increases, the ef- 
fective coupling of the superconductor to the Luttinger 
liquid renormalizes to zero. For a finite value of L and 
U > 0, the coupling is weak and one finds the usual 
Jisin</> Josephson relation. As L increases, J\ rapidly 
decreases and the transport of pairs through the chain 
vanishes in the L — > oo limit. 



III. THE EFFECTIVE SUPERFLUID WEIGHT 

£> 

In the following numerical study, we will be inter- 
ested in determining an effective superfluid weight D. If 
the pair phase varies linearly across a Hubbard chain of 
length L, then there will be a uniform Josephson current, 
and we will define D(L) by 



with <fio the phase difference across the Hubbard chain. 
The effective superfluid weight D is then given by D(L) 
as L — > oo. Here, we have set e = h = 1. The problem 
of determining D{L) is to create a linear phase change 
4>q/L across the Hubbard chain and then to measure j. 
The latter is straightforward since 

ji = — i[H,m] (6) 

so that for i ^ 1 or L, 

ji = -it^ ( c L c *+l^ - 4+l<r c i°) ■ ( 7 ) 
a 

At the boundary, when i = 1 (or L) we have to consider 
the boundary terms and add an extra current operator 

j'l = -iV(exp(icf))c\ 1 c\ l - exp(-j</>)c u ci T ) . (8) 

with a similar term for the right hand i = L boundary. 
The current density is independent of the position, and 
any of these expressions can be used with these end cor- 
rections to calculate j. 

The measurement of j is straightforward within the 
DMRG method^^. However, it is also necessary to es- 
tablish a uniform phase gradient. As noted in the previ- 
ous section, for a finite length L of the Hubbard chain, 
this can require tuning of the contact boundary pairing 
strength and the boundary scattering potential. Fortu- 
nately, for negative values of U, the contact interaction 
renormalizes to the perfect Andreev reflection fixed point 
as the length of the chain increases. However, when the 
finite system is doped away from half-filling, there are 
two parameters to tune and achieving a match such that 
the phase gradient over the length L is uniform becomes 
more difficult. For this reason, we have developed an ap- 
proach based upon extended contact interaction which 
will be discussed at the end of the next section. 

IV. RESULTS 

We use the Lanczos method for a system of size L = 8 
and DMRG for larger systems. The DMRG method is 
the standard finite-size algorithm, except for the use of 
complex numbers due to the arbitrary Josephson phases, 
and a special treatment of quantum numbers. The non- 
particle conserving boundary conditions mean that the 
total number of fermions cannot be used as a conserved 
quantum number. However, one can still utilize the num- 
ber of fermions modulo 2. This modulo-2 approach was 
first used inii. Within this approach the local pair-field 
A can take on a definite nonzero value. We have typically 
kept m = 200 states per block for the results presented, 
with a truncation error of about 10~ 8 . 

A. Point Contacts 



(5) 



In Figfl]we show Lanczos results for the Josephson cur- 
rent versus <j) — 4>r — 4>l through a half- filled Hubbard 
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FIG. 2: Effective superfluid weight D(L) of the half-filled 
Hubbard chain as a function of the Coulomb interaction U, 
for A = 1, V = 0, and chains of various lengths. We add for 
comparison the exact L — ► oo Bethe Ansatz results frorn^ for 
U< 0. 



<|)/Jt 

FIG. 1: Josephson current though a half-filled Hubbard chain 
with L — 8 as a function of the phase </> and for different 
values of the contact pairing A, and Coulomb interaction U, 
in units where the hopping t = 1. 

chain of L = 8 sites. For this half-filled, particle-hole 
symmetric case, with U < 0, the required site poten- 
tial Vi l = and the contact pairing strength A can 
be adjusted to achieve perfect Andreev reflection. For 
the non-interacting case, this is obtained for A/t = 1 as 
shown in the top panel of Fig[TJ For negative values of U, 
it is necessary to fine-tune A. When perfect Andreev re- 
flection is achieved, ji(4>) exhibits a sawtooth form with 
j(<fi) — D(L)cj)/L for — 7r < <fi < 7r. In this case, D(L) can 
be directly determined from j ((/>). For negative values 
of U , D rapidly approaches its asymptotic value when 
L 3> irt/\U\, so that the important requirement for de- 
termining D is to achieve perfect Andreev reflection at 
the ends. 

In Figure we show Lanczos and DMRG results for 
the superfluid weight D(L) of the half-filled chain for 
different values of the Coulomb interaction U. Here we 
have set V = and taken A = 1. The rcnormaliza- 
tion to perfect Andreev reflection is rapid for U < and 
the resulting effective superfluid weight D(L) varies little 
with L giving a value in close agreement with the exact 
Bethe- Ansatz results for the superfluid weight of the in- 
finite system, taken frorrA For U > 0, the system renor- 
malizes as L increases to the non-superconducting fixed 
point and the Josephson current is rapidly suppressed. 

In Fig[3]we show the DMRG results for the pair field 



amplitude along a chain at half-filling. The phase clearly 
varies linearly for negative U, while for positive values 
the modulus decays in a very short distance, a signature 
of the absence of superconductivity, and the phase order 
is disrupted by the small amount of noise in the DMRG. 

For the 8-site chain, we have seen that for the half- 
filled, particle-hole symmetric case it is necessary to tune 
the contact pairing strength A in order to achieve perfect 
Andreev reflection. For the non half-filled case, for finite 
L, there are two contact coupling parameters, A and V, 
that require tuning. 

In FigQ] we show DMRG results for the superfluid 
weight versus electron density and various values of the 
Coulomb interaction. Here, n is the electron density in 
the bulk of the chain, i.e. the center of the chain and 
far from the contacts). For comparison we show results 
for the superfluid weight D s for L — > oo obtained from 
Bethe Ansatz calculations^. For U = we have adjusted 
the values of A and V for maximum transmitivity, Eq. 2] 
For finite U we have set A = 1 and V\.l = 0, i.e. they 
are not optimized for perfect reflection. 

In FigfS] we show plots of D versus n for U = —2 
and chains of different lengths L. As L increases, D ap- 
proaches the exact result as the point contact bound- 
ary condition renormalizes to perfect Andreev reflection. 
However, to control this convergence it is in principle nec- 
essary to extrapolate the result to zero DMRG truncation 
error (large number of states m) and then take the infinite 
length limits. Hence, we would expect these curves to 
be more accurate if we were to fine-tune the parameters. 
However, this task has proven to be difficult. In order to 
overcome the difficulties of fine tuning the parameters in 
the Hamiltonian for optimal transmitance, we have stud- 
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FIG. 3: Pair field amplitude and phase along the L — 64 
chain for different values of U at half-filling. 
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FIG. 4: D as a function of the electron density for differ- 
ent values of the Coulomb interaction U. For U — 0, the 
boundary fields have been adjusted using Eqs.Q to achieve 
perfect Andreev reflection, while for U < we used A = 1, 
and V\,l = 0. We add for comparison the exact Bethe-Ansatz 
results from* (solid lines). 



ied the effects of using extended smooth contacts at the 
boundaries. 



B. Extended contacts 

In the previous section we have discussed a Hubbard 
chain of finite length L connected to superconductors 
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FIG. 5: Effective superfluid weigth D of a Hubbard chain 
with U — —2, connected to point contacts, as a function of 
density n, and for different lengths L. We add for comparison 
the Bethe Ansatz results in the thermodynamic limit (solid 
line) . 



through point contacts. We have seen that it is necessary 
to tune the pairing strength A and the boundary scatter- 
ing potential V in order to obtain a linear phase change 
along the chain. In this section, we explore the effects 
of extended contacts as an alternative way to eliminate 
the normal contact reflection. This technique is inspired 
by the smooth boundary conditions approach^. Here, 
we have applied the pair field end terms over a length i c 
on the end of each chain, with the coefficient dropping 
smoothly to zero as the distance from the end approaches 
l c . We have 



h.c. 



L 

A(L - i) (e^c^c^ + h. 



e=L- 



+ J2V(£) ni + V(L-e+l)n e (9) 

i=i e=L-e c +i 

where we take A(x) = A(l + cos(iE7r/4))/2 14 -. In the fol- 
lowing, we will set V{x) = and examine various widths 
£ c of the contact. 

In calculating the superfluid weight with the extended 
contacts one must utilize only the local properties in the 
center of the system. In particular, one must measure 
the current and the gradient of the phase in the center of 
the system. The phase varies linearly in the central re- 
gion of the chain, and this allows a numerical calculation 
of its gradient. It can also be shown that the effective 
superfluid weight can be extracted from the quantity 1 : 



J = 



j(x)dx = D{L)<t>. 



(10) 
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FIG. 6: Effective superfluid weight of a Hubbard chain (L — 
64) connected to smooth contacts of width £ c — 20, as a 
function of density n, and for different values of U. We add for 
comparison the Bethe Ansatz results in the thermodynamic 
limit (solid lines). 



In our calculations we simply replaced the integral by a 
sum over all the links. We find that the results obtained 
using the two approaches agree to within 1%. 

Figure [6] shows the results for the effective superfluid 
weight D for a Hubbard chain of length L — 64 with 
contacts of width £ c = 20. As in FigHJ the solid lines are 
the Bethe-Ansatz results for D s in the thermodynamic 
limit. As one can see, the DMRG results are in close 
agreement with the Bethe-ansatz results, except for n — 
1. It may be that logarithmic contributions affect the 
convergence of the DMRG for n = 1—. The extended 
contact approach provides a much closer match between 
the supeconducting leads and the Hubbard chain so that 
we have essentially achieved perfect Andreev boundary 



conditions without any need for tuning of parameters. 



V. CONCLUSIONS 

Here we have reported results of a numerical study of a 
one-dimensional Hubbard model coupled to external pair 
fields. DMRG calculations typically use open boundary 
conditions making it simple to couple the ends of an inter- 
acting system to a classical potential or magnetic fields. 
Here we have explored the numerics involved in coupling 
to a quantum pair field which can inject or remove pairs 
of electrons. We have seen how the pair transport varies 
as a function of the interaction U, the filling (n) and the 
length of the Hubbard chain. Various current-phase re- 
lations associated with the degree of Andreev reflection 
were clearly seen. A phenomenological effective super- 
fluid weight D was introduced and found to be in close 
agreement with Bethe-ansatz results for the superfluid 
weight of an infinite ring. 

Finally, the idea of an extended pair transfer contact 
was introduced. This was found to provide a useful way 
to effectively match the pair field injection such that the 
Andreev reflection approached unity. This is reminis- 
cent of the extended tapered connections used to match 
waveguides with different propagation characteristics and 
may prove useful in obtaining optimal matching of bulk 
leads to nanowires. 
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